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Abstract 


A robust and efficient computational method for reconstructing the 
three-dimensional displacement field of truss, beam, and frame 
structures, using measured surface-strain data, is presented. Known as 
“shape sensing”, this inverse problem has important implications for 
real-time actuation and control of smart structures, and for monitoring 
of structural integrity. The present formulation, based on the inverse 
Finite Element Method (iFEM), uses a least-squares variational 
principle involving strain measures of Timoshenko theory for stretching, 
torsion, bending, and transverse shear. Two inverse-frame finite 
elements are derived using interdependent interpolations whose interior 
degrees-of-freedom are condensed out at the element level. In addition, 
relationships between the order of kinematic-element interpolations and 
the number of required strain gauges are established. As cm example 
problem, a thin-walled, circular cross-section cantilevered beam 
subjected to harmonic excitations in the presence of structural damping 
is modeled using iFEM; where, to simulate strain-gauge values and to 
provide reference displacements, a high-fidelity MSC/NASTRAN shell 
finite element model is used. Examples of low and high-frequency 
dynamic motion are analyzed and the solution accuracy examined with 
respect to various levels of discretization and the number of strain 
gauges. 
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vector of the kinematic variables 
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vector of the approximated kinematic variables 

non vanishing components of the strain tensor for the Timoshenko 
beam theory in the ( x,y,z ) coordinate system 

components of the strain tensor for the Timoshenko beam theory 
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components of the stress tensor for the Timoshenko beam theory 
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experimentally evaluated section strains 
section forces and moments 

axial rigidity, shear rigidities, torsional rigidity and bending 
rigidities of the beam 

shear correction factors 

distributed loads along the x-, y-, and z-direction 
least-squares functional 

vector of the nodal degrees of freedom of the inverse finite 
element 

matrix of the shape functions relating the kinematic variables to 
the nodal degrees of freedom 

matrix of the shape functions relating the section strains to the 
nodal degrees of freedom 

weighting coefficients 

length, cross-sectional area and moments of inertia of frame 
inverse element 

number and axial coordinate of the locations where the section 
strains are evaluated 
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matrix and vector of the single inverse element governing equation 
matrix and vector of the whole structure governing equation 

linear Lagrange polynomials 

quadratic Lagrange polynomials 
special cubic polynomials 
quartic Lagrange polynomials 

orientation of a strain-gauge with respect to the beam axis 

external radius of the circular cross-sectio 

radius and thickness of the thin-walled circular cross section 

time varying applied force 

amplitude of the applied force 

frequency of the applied force 
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1 Introduction 


Real-time reconstruction of structural deformations using measured strain data, is a key technology for 
actuation and control of smart structures, as well as for Structural Health Monitoring (SHM) [1], Known 
as “shape sensing”, this inverse problem is commonly formulated with the assumption that multiple strain 
sensors at various structural locations provide real-time strain measurements. Most inverse algorithms use 
some type of Tikhonov’s regularization, which is manifested by constraint (regularity) terms that ensure a 
certain degree of solution smoothness (refer to [2-5] and references therein). 

Most of the shape sensing efforts focused exclusively on beam-bending problems. Davis et al. [6] used 
optimized trial functions and weights to reconstruct a simple static-beam response from discrete strain 
measurements. To model more complex deformations, the approach requires a large number of trial 
functions and strain sensors. Kang et al. [7] used vibration mode shapes to reconstruct a beam response 
due to dynamic excitation. In their approach, modal coordinates are computed using strain-displacement 
relationships and measured surface strains; the method requires the same number of mode shapes and 
strain sensors. Kim et al. [8] and Ko et al. [9] used classical beam equations to integrate discretely 
measured strains to determine the deflection of a beam. By regression of experimental strain data and by 
accounting for the applied loading, Kim et al. [8] obtained a continuous curvature function, leading to the 
evaluation of beam deflection. Ko et al. [9] developed a load-independent method by approximating the 
beam curvature using piece-wise continuous polynomials; the authors demonstrated the validity of a one- 
dimensional scheme by evaluating the deflection and cross-section twist of an aircraft wing. 

To enable shape-sensing analyses of plates undergoing bending deformations, Bogert et al. [10] 
examined a modal transformation method that allows the development of suitable strain-displacement 
transformations. The approach makes use of a large number of natural vibration modes. When applied to 
high-fidelity finite element models, however, the method requires a computationally intensive eigenvalue 
analysis and a detailed description of the elastic and inertial material properties. Jones et al. [11] 
employed a least-squares formulation for shape sensing of a cantilever plate, where the axial strain was 
fitted with a cubic polynomial. The strain field was then integrated with the use of approximate boundary 
conditions at the clamped end to obtain plate deflections according to classical bending assumptions. 
Shkarayev et al. [12,13] used a two-step solution procedure: the first step involves the structural analysis 
of a plate/shell finite element model, and the second, a least-squares algorithm. The methodology 
reconstructs the applied loading first, which then leads to the displacements. In a series of four papers, 
Main 5 on and co-authors [14-17] developed a finite element formulation that seeks the solution for the 
displacements and loads simultaneously, requiring a priori knowledge of a subset of applied loading and 
the material properties. The solution procedure minimizes a cost function consisting of unknown loads 
and differences between the measured and estimated quantities (displacements or strains); the cost 
function is regularized by way of equilibrium constraints. The number of unknowns is three times the 
number of the degrees-of-freedom in the finite element discretization. Importantly, the accuracy of the 
solution strongly depends on the choice of suitable weights; these are computed from a complex 
procedure involving the probability distributions of the unknown loads and measured data. In [16,17], 
sensitivity analyses were carried out for truss structures, investigating variations in the input data as well 
as the modeling errors. Nishio et al. [18] employed a weighted-least-squares formulation to reconstruct, 
on the basis of measured strain data, the deflection of a composite cantilevered plate. The weighting 
coefficients in the least-square terms were adjusted in order to account for the inherent errors in the 
measured strain data. The weights were computed for a given data-acquisition apparatus, load case, and 
test article, with the consequent difficulties in generalizing the procedure. 

Many of the aforementioned inverse methods either lack generality with respect to structural topology 
and boundary conditions, or require sufficiently accurate loading and/or elastic -inertial material 
information - the kind of data that are either unavailable or difficult to obtain outside the laboratory 
environment; for these reasons, such approaches are generally unsuited for use in onboard SHM 
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algorithms. A well-suited algorithm for SHM should be: (1) general enough to accommodate complex 
structural topologies and boundary conditions (e.g., built-up aircraft structures), (2) robust, stable, and 
accurate under a wide range of loadings, material systems, inertial/damping characteristics, and inherent 
errors in the strain measurements, and (3) sufficiently fast for real-time applications. 

An algorithm that appears to fulfill the aforementioned requirements was recently developed by 
Tessler and Spangler [1,19]. The methodology, labeled the inverse Finite Element Method (iFEM), 
employs a weighted-least-square variational principle which is discretized by CO-continuous finite 
elements that accommodate arbitrarily positioned and oriented strain-sensor data. The iFEM framework, 
providing accurate and stable solutions of the displacement and strain fields for a discretized structural 
domain, is amenable to any type of structural modeling including frame (truss and beam), plate, shell, and 
solid idealizations. Because only strain-displacement relations are used in the formulation, both static and 
dynamic response can be reconstructed without any a priori knowledge of material, inertial, loading, or 
damping structural properties. To model arbitrary plate and shell structures, Tessler [20] developed, using 
first-order shear-deformation theory, a three-node inverse shell element. Numerically generated [20] and 
experimentally measured-strain data [21,22] were used to assess robustness and accuracy of the 
formulation. 

The present paper consolidates the authors’ recent efforts [23-24] and presents the development and 
assessment of simple and efficient inverse-frame finite elements. The methodology permits effective and 
computationally efficient shape-sensing analyses to be performed on truss, beam, and frame structures 
instrumented with strain gauges. The kinematic assumptions are those of Timoshenko shear-deformation 
theory [25]; they incorporate stretching, torsion, bending, and transverse shear deformation modes in 
three dimensions. The formulation uses a least-squares variational principle that is specialized from [19] 
for three-dimensional frame analysis. The variational framework, in conjunction with suitable finite 
element discretizations involving inverse finite elements, yields a system of linear algebraic equations; the 
equations are efficiently solved for the unknown displacement degrees-of-ffeedom (dof s), thus providing 
the deformed structural-shape predictions. 

In the remainder of the paper, the kinematic assumptions for a three-dimensional frame are discussed, 
followed by the description of the least-squares variational principle suitable for three-dimensional 
deformations of frame structures. This is followed by a discussion of two C°-continuous, inverse-frame 
elements that use the well-established interdependent interpolations that resolve the shear locking effect. 
Finally, to examine the predictive capabilities of the inverse elements for a given set of distributed strain 
gauges, shape-sensing studies are carried out for a cantilevered beam undergoing harmonic excitations in 
the presence of structural damping. 

2 Governing equations 

Consider a straight frame member of constant cross-section referred to the three-dimensional 
Cartesian coordinates ( x,y,z ) as depicted in Figure 1; the coordinate origin, O, is located at the cross- 
section’s center of mass, which is also coincident with the shear center. The longitudinal, elastic x-axis is 
normal to the cross-sectional plane (y, z), where y and z are the cross-section’s principal inertial axes. The 
frame member has length L and its cross-section has area A, area moments of inertia with respect to the 
y - and z -axis 7 v and I _ , respectively, and polar moment of inertia I p = 7 v + 1 _ (Figure 1). The frame 

member is made of an isotropic homogeneous material, represented by the elastic constants: E (Young’s 
modulus), G (shear modulus), and v (Poisson ratio). 

The three Cartesian components of the displacement vector that are consistent with the kinematic 
assumptions of Timoshenko theory [25] for three-dimensional deformations are given by 
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( 1 ) 


u x (x, y, z ) = u (x) + z6 y (x) - yO, (x) 
u y (x, y, z) = v(x) - z6 x (x) 
u z (x, y, z) = w(x) + yO x (x) 


where u x , u y , and u are the displacements along the x , y , and z axes, respectively, with u , v , and 

w denoting the displacements at y = z = 0 ; 0 X , 9, . and 0 are the rotations about the three coordinate 

axes; positive orientations for the displacements and rotations are depicted in Figure 1 . These kinematic 
assumptions neglect the effect of axial warping due to torsion, i.e., each cross-section remains flat and 
rigid with respect to thickness-stretch deformations along the y and z axes. The six kinematic variables 
can be grouped in vector form as 


U = { M , v, w, 0 X , 0 y , 0^ T (2) 



Figure 1 : Beam geometry and kinematic variables. 


Adhering to the small-strain hypothesis, the non-vanishing strain components have the form 

£ x (x, y, z) = e, (x) + ze 2 (x) + ye 3 (x) 
r xz (x,y) = e 4 (x) + ye 6 (x) (3) 

yjx, z) = e 5 (x)-ze 6 (x) 


where 


6 


(4) 


e(u) = {< 


P 


" 2 ’ 


- 3 ’ 


"4’ 


"5 ’ 




denote the section strains of the theory, given by 


e x (x) = u x (x) e 4 (x) = w x ( x ) + 6 y (x) 

e 2 (x ) = O y Jx) e 5 (x) = v x (x)-O z (x) (5) 

e 3 (x) = -0 zx (x) e 6 (x) = 0 xx (x) 


The section forces ( N , Q y , and Q z ) and moments (M x , M y , and M.) are related to the section strains, e t , by 
way of the constitutive equations (refer to Figure 2) 


N = A A M , = J , e 6 

Qv = G y e s M v = D y e 2 

Q = G z e A M = D z e 3 


(6) 


where A x = EA is the axial rigidity; G y = k y GA and G = k_GA are the shear rigidities, with k y and 
kz denoting the shear correction factors; J x = GI P is the torsional rigidity and D = El and D = El _ 
denote the bending rigidities. 


M, 


1 



Figure 2: Beam section forces and moments. 


The equilibrium equations that correspond to the distributed loads q t (x) , q X (x) , and q (x) along the x, 
y and z directions, are 


7 


dM 


dN _ 

dx; 

dQ^ = 

dx; 

dQ^ = 

dx; 


a x 

% 

Vz 


dx; 

dM 


= 0 


dx; 

dM 


dx; 


L = Q, 

= Q, 


( 7 ) 


To reconstruct the deformed shape of a frame-member for which certain in-situ strain measurements 
are known, a functional ®(u) that matches in a least-squares sense the complete set of analytic section 

strains, e(u) , to the in-situ section strains, e E , is minimized with respect to the kinematic variables, u ; 
®(u) functional can be written as 


®(u) = e(u)-e £ 


( 8 ) 


This functional may be used in a finite element framework by introducing a discretization in which the 
element kinematic field is interpolated by C° -continuous shape functions, 


u(v) - u* = N(.x;)iT 


(9) 


where N(x) denotes the shape functions and u e nodal dof s. Thus, the total least-squares functional is a 

N 

sum of the individual element contributions, ®‘ (u / ) , i.e., ® = ^ d>' , with N denoting the total number 

e=l 

of elements. Accounting for the axial stretching, bending, twisting, and transverse shearing, the element 
functional is given by the dot product of the weighting coefficient vector, 

w = {w jt } = {w[ ) , w° 2 (l e y /A e ^, wl(rjA e ), w”, , and the least-squares 

component vector, ® = |®^ | , 

d>V) = w-<I> (10) 


where (k = 1, 6) denote dimensionless weighting coefficients; A e , V , /f, and I e p are, 

respectively, the cross-section area, moments of inertia with respect to the y - and z -axis, and polar 
moment of inertia of the element, and 
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2 


(k = 1, 6) 


( 11 ) 


O 


e 

k 


11 


n 


X[ g * (*<■)-<'] 


denote the least-squares components of the element functional, where L e denotes the element length; n 
and x i ( 0 < x i < 111 ) are, respectively, the number and the axial coordinate of the locations where the 
section strains are evaluated, and the superscript si is used to denote the section strains that are computed 
from the strain-sensor values (experimental values) at the location x i . The w° k coefficients may be 
assigned different values to enforce a stronger or weaker correlation between the measured section-strain 
components and their analytic counterparts, i.e., a larger value of w]' enforces a stronger correlation, 
whereas a smaller value enforces a weaker correlation. 

Substituting Eq. (9) into Eq. (3) gives the section strains in terms of the nodal dof s as 


e k =B k u e 6 ) 


( 12 ) 


A vector form of Eq. (12) that incorporates all six section strains is given by 


e(u) = B(x)iT 


(13) 


where the matrix B ( xj contains the derivatives of the shape functions N(x) ( B ( xj is defined in 

Appendix B for the case of shape functions presented in Sec. 3.) Substituting Eq. (13) into Eq. (11) and 
then Eq. (10) results in the following quadratic form 


<D e =|(u f fkV-(u e ) T f'+c' 


(14) 


where c e is a constant while k and f are defined as follows 


k e 




(15) 


with 
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(k = 1, 6) 


(16) 


k 


e 

k 


L e 


n 


X[b[N)b ( (*,)], f; = 


L e 


n 


Z[b[(*,K] 


Note that k resembles an element stiffness matrix of the direct finite element method and f resembles 
the load vector; k e is a function of the measurement locations, x ; , whereas f possesses the measured 

strain values. Minimization of functional d>‘ (see Eq. (14)) with respect to u leads to the inverse 
element matrix equation 


kv=r 


(17) 


The assembly of the finite element contributions, while accounting for the appropriate coordinate 
transformations and by specifying problem-dependent displacement boundary conditions, results in a non 
singular system of algebraic equations of the form 


KU = F 


(18) 


The solution of these equations for the unknown dof s is efficient: the K matrix is inverted only once, 
since it is independent of the values of the measured strains. The F vector, however, is dependent on the 
measured strain values that change during deformation. Thus, at any strain-measurement update during 
deformation, the matrix-vector multiplication provides the solution for the unknown nodal displacement 
dof s, U = K 1 F, where K 1 remains unchanged for a given distribution of strain sensors. 

The remaining part of the element formulation involves the selection of suitable shape functions, 
symbolically defined by Eq. (9), and the computation of the experimental section strains, e “ , appearing 
in Eqs. (12), (13). In Sec. 3, the shape functions for two alternative inverse-frame elements, each having 
two nodes and twelve dof s, are derived. In Sec. 4, a procedure for computing e“ is described; it relates 
the number of strain gauges to the interpolation order of the shape functions. 

3 Element shape functions 

In this section, inverse frame elements of 0 th - and l st -order are formulated. The elements use C°- 
continuous interdependent interpolations that enable excellent predictions even for very slender frame 
members, without incurring any form of excessive stiffening due to shear locking [26]. The 0 lh -order 
shape functions are guided by Timoshenko equilibrium equations, Eqs. (7), that correspond to the forces 
and moments applied exclusively at the end nodes, resulting in constant distributions of the transverse- 
shear section strains. The l st -order shape functions accommodate Eqs. (7) for uniformly distributed 
transverse loads, giving rise to linear distributions of the transverse-shear section strains. 

A frame element is referred to a local axial coordinate xe[^0, Z) J, where Zf denotes the element 

length. Furthermore, a non-dimensional coordinate £, = (2x1 L — 1) e [-u] is used to define the 
element shape functions (Figure 3). The initial nodal configurations are defined by the two end nodes, 1 
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(at £ = — 1) and 2 (at £ = +1) and one or three interior nodes. Thus, the initial configuration for the 0 th - 
order element has the interior node, r (at the midspan, % = 0 ); whereas the interior nodes of the l st -order 
element are <7 (at E, = — 1 / 2 ), r (at % = 0 ), and 5 (at E, = + 1 / 2 ) . 



1 q r s 2 


/ 

— c^— 

—^3 

L e 

—O— 

\ 





/ 


Figure 3: Inverse finite element geometry and nodal topology. 


The initial nodal configurations of the 0 th - and F'-order elements are readily reduced to two nodes and 
twelve dof s by condensing out the interior dof s at the element level in a manner analogous to static 
condensation. The resulting elements have three-displacement and three -rotation dof s at each end node 
(Figure 4); thus, the dof s vector of the elements is 


U — jwp t ; j, Wp 0 y \ , 0y\i 0 7 \ t ^2’ @x2 ’ ^v2’ ^ 2 } 


u 


u 

V 


V 

w 


w 

o x 


o x 



0 y 

kJ 

1 

kJ 


Figure 4: Two-node inverse finite element. 


(19) 


The process of condensing out the interior dof s results in the reduced element equations, k' A ,u/ = f/ , 

where k' A , is a function of the partitioned parts of the original k matrix, and u/, contains the end-node 

dof s. Since the unreduced k matrix is independent of the strain values, so is the k' A . matrix. This 
implies that even for the elements with the condensed-out interior dof s, the corresponding system matrix, 
K , is strain-value independent (refer to Eq. (18).) 

3.1 0 th -order element 

The formulation of the 0 th -order element is guided by Eqs. (7) for the loading case of end-node forces 
and moments. For this case, the axial force, twisting moment, and shear forces are constant along the 
element; whereas the bending moments are linear. Eqs. (7) in terms of the section strains (after Eqs. ( 6 ) 
have been introduced) indicate that the section strains e, (i= 1,4-6) are constant, and e, (i=2,3) are linear. 
From Eqs. (5), it is deduced that u and 0 , are linear, 0 and 0 are parabolic, and v and w are cubic. 
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The inter-relationship of the polynomial order of the deflection variables, v and w , and bending 
rotations, 0 and 0 , can also be inferred from the definitions of the transverse shear section strains [24] 

e 4 =[w x +d y ), e 5 =(v x -0 7 ) (20) 

It is understood that in e 4 both w v and 6 should be represented by the same polynomial order and, 
similarly, in e 5 both v t and 6 should also have matching polynomial orders. The above interpolations 

give rise to quadratic interpolations for e A and e 5 ; they permit a consistent reduction of interior dof s for 
v and w by requiring a constant variation of these section strains across the element span 

e 4 = const., e 5 = const. (21) 

The complete set of interpolations for this element is thus given by 


u {^)=H L ?{^) u o 

1 = 1,2 

*,(£)= I L ?(mr 


0d£)=Z4’(£R 

*= 1,2 

<>M)= z 


(22) 


j= Ur, 2 


j=l,r,2 


’(f) = Z4"(4K- Z «'(f)=Z4 > (f)’*' l + Z N T(?)0„ 


i= 1,2 


7=1, t ,2 


i=l,2 


7=1, r ,2 


where Lf } (<^) (i = l,2) and Z)/' (<^) (j =l,r,2) are, respectively, linear and quadratic Lagrange 
polynomials, and N { ^ (<g j (7 = 1, r , 2) are special-form cubic polynomials (refer to Appendix A.) Static 
condensation can be used to condense out the two interior dof s (0 and 6 _ r , refer to Eqs. (22)), thus, 
achieving a two-node element with twelve dof s (Figure 4, Eq. (19)). 

3.2 l st -order element 

Consider a frame element loaded by uniformly distributed transverse loads, q (x) and q. (x) . From 
Eqs. (7), after the substitution of Eqs. (6), it is readily deduced that e 4 and e. need to be linear and e, 
and e 3 parabolic. The a and 6 X variables remain linear as in Eq. (22). Moreover, v and w are quartic 
whereas 0 y and 6 are cubic. Following the constraint strategy for e 4 and e 5 , these section strains are 
set to be linear for this element. The resulting interpolation polynomials are given by 
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(23) 


'(£) = ZC(£K 


4 60 = Z 4" (<?K 


i=l,2 


i=l,2 


*■(£)= Z 4 4, («)v,„ 

£=!,#, r,s, 2 


W 


60 = Z C(fl 

fc=l,#,r,s,2 


w. 


460 = Z4”604< + Z *“(£K *,(£)= Z4‘’(fR- Z *"'(£) 


(3) 


/= 1,2 


£=l,g,r,.s,2 


i=l,2 


k=l,q,r,s,2 


— (3) 

where Nk (c ) = 1, q, r, s, 2) are cubic polynomials that satisfy the conditions 


e 4 = linear , e 5 = linear (24) 

For the detailed expressions of the Nk (c ) polynomials, refer to Appendix A. The interior dof s (v , 
v r , i’ t , w , w r , and w s , refer to Eqs. (23)) are condensed out at the element level, leading again to a 
twelve dof s inverse element as in Eq. (19). 

4 Input data from surface strain measurements 

A key step in the iFEM formulation is to evaluate the section strains due to experimentally measured 
surface strains. In this section, the relationships between the measured surface strains and the six section 
strains, e h are established. Also discussed are the strain gauge positions along the frame axis and their 
angular orientations that enable the complete description of the experimental section strains. For 
illustration, the present analysis is restricted to frame members with circular cross-sections; the adopted 

cylindrical coordinate system (0, x, r) is shown in Figure 5. 



Figure 5: Orthogonal and cylindrical coordinate systems. 

4.1 Section strains derived from linear strain gauge measurements 

Taking the usual assumption of negligible a v and a _ , then cr t and z xd are the only non-zero stress 

components acting on the external surface r = R ext (Figure 6(a)). The corresponding strain state, 
represented in Figure 6(b), is 
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(a) Stress state. 


(b) Strain state. 


Figure 6: Stress and strain states on the frame external surface (t-R ex t) in the cylindrical 

coordinate system. 

Consider a linear strain gauge placed on the external surface at x = x ; , at a particular 6 and with an 
angle (3 with respect to the beam axis (Figure 7); is a local Cartesian coordinate system 

having its origin on the frame external surface, the x* -axis along the strain gauge measurement axis ( s * 2 ), 
the x* -axis on the frame surface and the x 2 -axis normal to the frame surface and coincident with r-axis. 



Figure 7: Location and coordinate system of a linear strain gauge placed on the frame 

external surface. 

Using appropriate strain-tensor transformations from the (0, x, r) to [x*, xj , x! j coordinates [27], the 
• • . * 

relationship between the measured strain s 7 and the strain tensor components in Eqs. (25) becomes 


s* 2 = s x cos^ f3 +s 0 sin 1 p + y x0 cos p sin p 
or, using the second of Eqs. (25), 


(26) 
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s* 2 = £ x (cos 2 P - v sin 2 /?) + y x9 cos p sin p 


(27) 


Substituting r = R ext in Eqs. (3), yields 

= e i + e 2 R », sin 0 + e 3 R exl cos 0 

(Z Oj 

r x e = e 4 cos 0 - e 5 sin 0 + e ( R ext 

. * • 

Substituting Eq. (28) into Eq. (27) results in the relation between the measured strain s 2 and the six 
section strains at x = x t 


s* 2 (x i ,0,p) = e i (x,)(c 2 -V5 2 ) 

+e 2 (x i )(c 2 p -vs 2 fi )s 0 R ext 

+e 3 {x i )(c 1 lj -vs 1 p )c 0 R ext (29) 

+e 4 {x i )Cf j s p c g 

-e 5 (x t )c p s p s g 

+e 6 (x i )c fr s /j R ext 

where c e = cos 0 . s g = sin 0 , c„ = cos p , and s f = sin p . 

4.2 Strain gauge distributions 

The iFEM formulation minimizes, in the least-squares sense, Eq. (8), where e E are the section strains 
computed from the measured strains. Thus, an important question arises what constitutes the minimum 
number of strain measurements. 

For the O th -order element, e 3 , e 4 , e 5 and e (i are constant, whereas e 2 and e 3 are linear with respect to the 
axial coordinate, x; these section strains necessitate eight strain measurements. Similarly, for the l st -order 
element, e x and e 6 are constant, e 4 and e 5 are linear, and e 2 and e 4 are quadratic, thus requiring twelve 
strain measurements. 

A further reduction of strain measurements is possible if one invokes Eqs. (7). Substituting Eqs. (6) 
into Eqs. (7) results in 


D v e 2, 


= G z e 4 , 


D e 4 

z 3,x 


= G y e s 


(30) 


For the 0 th -order element, Eqs. (30) give rise to two constraints equations, thus reducing the minimum 
number of strain measurements to six; whereas, for the l st -order element, Eqs. (30) give rise to four 
constraint equations, thus reducing the minimum number of strain measurements to eight. It is worth 
noting that this procedure should be viewed as a convenient means of reducing the required number of 
strain gauges by solving for e 4 and e 5 analytically rather than measuring these quantities experimentally. 

Since the strain gauges can be placed anywhere along the beam surface, the distributions considered 
in this study are summarized in Table 1 (also refer to Figures 8 and 9). To refer to a specific combination 
of the element type and strain gauge configuration, a compact notation, #-#E, is used; where the first 
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position, #, refers to the element order (0 or 1), the second position, #, indicates the number of strain 
gauges per element (6 or 8), and the letter “E” indicates that Eqs. (30) have been used in the formulation. 
The strain gauges are placed at different positions x = (L73, L72, 2L73) along the element. The strain 
gauge angular orientations (0,j3) are also allowed to be different; for example, ( 0, [i ) = (-2rc/3, tc/ 4) 
indicates that the strain gauge is placed at the circumferential angle 0=-2n/3 and is oriented with an angle 
/?=7t/4 with respect to the frame x-axis (Figure 7). 

Table 1 : Strain gauge distributions x, ( 0, (0 ) corresponding to the 0 th - and 1 st -order elements. 


Element- 
strain gauge 
notation 

Orientation (0,/3) of 
strain gauges at 
x=L73 

Orientation (0,j3) of 
strain gauges at x=L e H 

Orientation (0,j3) of 
strain gauges at 
x=2L73 

0-6E 

- 

(-271/3,0), (-271/3,71/4), 
(0,0), (0,7t/4), 
(271/3,0), (27t/3,7i/4) 

- 

1-8E 

(-271/3,71/4) 

(-2tx/ 3,0), (-27i/3,7t/4), 
(0,0), (0,71/4), 
(271/3,0), (27i/3,7i/4) 

(27i/3,7i/4) 



Figure 8: 0-6E strain gauge distribution. 



Figure 9: 1-8E strain gauge distribution. 
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5 Numerical results 


A simple cantilevered beam subjected to dynamic loading is analyzed to assess the accuracy of the 
inverse finite element method. The beam is made of an aluminum alloy (£=73,000 MPa, V =0.3, and 
p=2557 Kg/m 3 ) and has a thin-walled circular cross-section with the average radius R =39 mm, wall 
thickness s = 2 mm, and length L=800 mm. The beam is fully clamped on one end and subjected at the 
other end to a harmonic vertical force F z (t) (where t denotes time) acting in the z-direction at frequency f 0 , 
i.e., 


F z (t) = F z o sin (2 nfy) 


(31) 


where F_ () is the force amplitude (F 0 = 10 3 N.) To simulate the experimental-strain measurements and to 

assess the accuracy of the inverse method, high-fidelity direct FE analyses are performed using QUAD4 
shell elements in MSC/NASTRAN. The model consists of 1 14 elements along the cross-sectional 
circumference and 360 elements along the beam axis, for a total of 41,040 elements and 41,156 nodes. 
The tip force is applied at the cross-sectional center at a node which is connected to all other nodes within 
the cross-section by means of multi -point constraints (or MPC’s) [28], 

The dynamic response of the beam is calculated using a modal transient analysis in MSC/NASTRAN 
and keeping the modes up to 5,000 Hz, with the inclusion of viscous damping of magnitude 5% with 
respect to the critical value at each frequency. In the frequency range from 0 to 5,000 Hz, 51 modes are 
present: these include the first lowest flexural beam modes, 1F-5F, appearing twice due to the cross- 
section symmetry, and the first membrane mode (1M). Table 2 summarizes the order of the global modes, 
their type, and corresponding frequency value/; the first three flexural mode shapes are shown in Figures 
10-12. The other modes in this frequency range are associated with shell modes describing cross-sectional 
distortion and are neither shown nor tabulated. 


Table 2: Global modes of the cantilevered beam in the frequency range of 0-5,000 Hz. (F-type modes are 

flexural; M-type modes are membrane.) 


Mode order 

1 st and 2 nd 
modes 

3 rd and 4 th 
modes 

12 th mode 

13 th and 
14 th modes 

30 th and 
31 st modes 

40 th and 
41 st modes 

Mode type 

IF 

2F 

1M 

3F 

4F 

5F 

/ frequency 
[Hz] 

126.8 

729.5 

1,670 

1,835 

3,187 

4,671 
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Figure 10: 1 st flexural mode (IF, / = 126.8 Hz). 



Figure 1 1 : 2 nd flexural mode (2F, / =729.5 Hz). 



Figure 12: 3 ld flexural mode (3F, / =1,835 Hz). 


To investigate the accuracy of the iFEM modeling for dynamic applications in both low- and high- 
frequency regimes, three different values of the applied-force frequency f 0 have been considered, namely: 
fo = 60 Hz (about half of the fundamental frequency), / 0 =450 Hz (halfway between the IF and 2F modes), 
and f 0 = 1,400 Hz (halfway between the 2F and 3F modes). Figures 13-15 compare the tip-deflection time 
histories, w max (t), calculated by means of the high-fidelity FEM shell model using MSC/NASTRAN and 
the corresponding iFEM frame-element models. The tip deflection of the NASTRAN model corresponds 
to the cross-sectional center, and is computed at a node which is connected to all other nodes within the 
cross-section by means of MPC’s. The present iFEM models used the strain-gauge distributions in Table 
1 and the uniform weight coefficients w° k = 1 (k = I, ..., 6) in Eq. (10); the strain values were taken from 

the nodes (at the specific locations in Table 1) of the NASTRAN model. It is noted that slightly more 
accurate strain values reside at the element Gauss points. However, considering the high fidelity of the 
reference FEM model, the “measured” strains taken at the nodes are quite satisfactory. 

For the low-frequency loading of / 0 =60 Hz, a single 0 th -order inverse element gives accurate results, 
with a maximum error in the tip deflection of 2.3% (Figure 13). At this excitation frequency, when 
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t > 0. 1 , viscous damping has reduced the structural vibrations to a steady state response, proceeding at a 
constant amplitude and the same frequency as the forcing function. When the excitation frequency of the 
forcing function is increased, the response has a longer transient region, which is manifested by 
interactions between the natural modes of vibration and that due to the applied dynamic loading. To 
model the transient response at higher frequencies, finer discretizations are required. Thus, for f 0 = 450 
Hz, a two-element, l st -order model yields a 1.1% error in the maximum deflection (Figure 14). At the 
fo- 1,400 Hz frequency, a three -element iFEM discretization using the l st -order element results in the 
maximum deflection error of 2.0% (Figure 15). These results clearly demonstrate that the methodology is 
highly efficient, requiring only few inverse elements and strain gauge measurements, and is applicable not 
only for the steady state portion of the response but also for the transient regime at high frequencies. 



0 0.05 0.1 0.15 

t (sec) 


Figure 13: Tip deflection w max of the beam loaded by a transverse concentrated force F z at 

fo= 60 Hz. 
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Figure 14: Tip deflection w max of the beam loaded by a transverse concentrated force F z at 

fo = 450 Hz. 



Figure 15: Tip deflection w max of the beam loaded by a transverse concentrated force F z at 

fo= 1,400 Hz. 

2( 



6 Conclusions 

In search of a suitable computational method for use in Structural Health Monitoring (SHM) systems, 
an inverse Finite Element Method (iFEM) has been formulated to perform the displacement- 
reconstruction analysis (shape-sensing) of three-dimensional frame structures undergoing static and 
dynamic deformations. The methodology uses a least-squares variational principle, which is discretized 
by CO-continuous displacement-based inverse frame elements. Linear strain-displacement relations and 
their components, known as section strains, are based on the Timoshenko (first-order) shear deformation 
theory that includes the deformations due to stretching, torsion, bending, and transverse shear. The 
variational statement enforces experimentally measured strains to be least-square compatible with those 
interpolated within the inverse frame elements. The implementation of this least-square compatibility is 
accomplished using the individual section strains. 

Two inverse frame elements, each having two nodes and twelve dof s, have been developed. The 0 th - 
order element has a constant shear-section strain along the element length, whereas the l st -order element 
has a linear shear section strain. The element shape functions are based on interdependent interpolations 
that ensure locking-free bending of slender frame members. The element interpolation order is linked to 
the definition of the number and orientation of the uniaxial strain gauges that are necessary for the 
analysis. Two simple and effective strain gauge distributions have been selected and used in the numerical 
examples. 

The present shape-sensing capability has been demonstrated on a thin-walled, circular cross-section 
cantilevered beam subjected to harmonic excitations in the presence of structural damping. To provide the 
simulated strain-gauge measurements, as well as the reference displacements, a high-fidelity shell finite 
element model was developed using the MSC/NASTRAN commercial code. Low- and high-frequency 
dynamic beam motions were analyzed and time history of the tip deflection examined, comparing several 
iFEM discretizations and strain-gauge schemes. The iFEM shape-sensing analysis, which is based only on 
the strain-displacement relations and the measured strain data (without any reliance on the material, 
inertial, or damping properties of the structure), has been shown to be highly effective and efficient in 
predicting the dynamic structural response of a damped beam. Accurate predictions of both the steady- 
state and transient response required only few elements and strain-gauge measurements, where the higher- 
frequency excitations necessitated somewhat higher fidelity of the iFEM models. 

Although beyond the scope of the present effort, additional studies need to be performed, including: 
(a) shape-sensing analysis of spatial frame structures using the strains measured in a laboratory, and (b) 
studies of the strain-gauge distributions that provide optimal (or nearly optimal) solutions. 
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Appendix A 

The 1 st , 2 nd , and 4 th -degree Lagrange shape functions are given as 


• 1 st degree 

(ai) 

• 2 nd degree 


• 4 th degree 


(A2) 


[C,4*>] = D( 4 f- 1)[(#-1),(£ + 1)] 

(A3) 

[i« > ,Lr.4 4l ] = ^(l-f)[4f(2f-l),3(l-4f),4f(2f + l) 


where E, = 2x1 Zf — 1 e [— 1 , 1 ] is a non-dimensional axial coordinate; x e j^0, L J and If denotes the 

element length. The subscripts 1 and 2 represent the end nodes, whereas q, r, and s denote the uniformly 
spaced interior nodes. 

The 3 rd -degree shape functions, TV' 3 ’ (fj , of the 0 th -order element have the form 
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<A4) 


— ( 3 ) 

whereas the cubic Nk (c) shape functions of the l st -order element are 
N? ^ M\n? + + (A5) 

- —1 Jx_/ 

Appendix B 

The matrix B relating the section strains to the element dof s may be written as follows 

B, B n 0 B l2 

B 2 b 21 b 2c b 22 

b 3 _ B 31 b 3c b 32 

b 4 b 41 b 4c b 42 

B 3 B 51 b 5c b 52 

®6_ b 61 0 b 62 
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b*=^[ 43 o] 


B* = 2[o -41 


B 4o = 


4i’+- 


EE 


(2) A 


11 I 


0 - 


+ 


EE 


(2)1 


(B3) 


For a l st -order element, 0 is a 1x6 null matrix and the other sub-matrices are defined below 


B „=fr[4] 0 0 0 0 0] 


B, = 


TT(3) 


0 0 0 I™ 0 


B 3; s — 
3 E 


0 Nfl 0 0 0 -ig 


B 4i - 


B 5i - 


0 0 


l {4) + 
1,% + 


L‘N E 


E’L (2) 

0 i- 0 


l (4) + 


L e N 


(3) 1 


0 ( 2 ) 


0 0 0 - 


EE 


B,^-[0 0 0 Lfl 0 0] 


0‘ = 1,2) 


(B4) 


b 2c - 


n 77 (3) n T7 (3) o 77 (3) 

0 N q,t; 0 N r,4 0 Ns,4 


B, = 

3c L e L 
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' — (3) — (3) — (3) 

N Chi 0 N r,s 0 Ns,e 0 


B 4c - 
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